Density matrix renormalization group for disordered bosons in one dimension 
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We calculate the zero-temperature phase diagram of the disordered Bose-Hubbard model in one 
dimension using the density matrix renormalization group. For integer filling the Mott insulator 
is always separated from the superfluid by a Bose glass phase. There is a reentrance of the Bose 
glass both as a function of the repulsive interaction and of disorder. At half-filling where no Mott 
insulator exists, the superfluid density has a maximum where the kinetic and repulsive energies are 
about the same. Superfluidity is suppressed both for small and very strong repulsion but is always 
monotonic in disorder. 



The interplay between disorder-induced localization 
and interactions has attracted a great deal of attention in 
recent years. The simplest model including both aspects 
in a nutshell is a Hubbard model with random site en- 
ergies and a local repulsive interaction for either bosons 
or fermions with opposite spin. Unfortunately there are 
essentially no analytical results for this model if both 
disorder and interactions are present, not even in one di- 
mension. For Id bosons, however, there is a weak disor- 
der, perturbative calculation by Giamarchi and Schulztl, 
who found that the superfluid ground state with quasi 
long range order survives disorder up to a critical point, 
where the effective exponent K in the decay of the one 
particle density matrix is equal to 2/3. More generally, 
the qualitative physics of the Bose- Hubbard model in any 
dimension, and in particular the scaling behaviour near 
critical points has been elucidated by Fisher et al.cl. For 
quantitative results, however, it is necessary to resort 
to numerical simulationstrcl. The latter were performed 
using path integral (or "world line") Monte Carlo cal- 
culations which become increasingly difficult in the most 
interesting limit of zero temperature. Now at least in one 
dimension there is an inherently zero temperature numer- 
ical technique for interacting quantum problems, the den- 
sity matrix renormalization group (DMRG) method de- 
veloped by WhiteO. It is therefore natural to try employ- 
ing this method to the disordered Bose-Hubbard model 
in one dimension. This was first done by Pai et al.u, 
who calculated the associated phase diagram for integer 
filling. As expected, it exhibits a Mott insulating, a su- 
perfluid and also a Bose glass phase, the latter appearing 
only for sufficiently strong disorder. Quite recently, how- 
ever, their results were seriously questioned by Prokof 'ev 
and Svistunov, who performed rather precise quantum 
Monte Carlo calculations!! Based on that, it was argued 
that the DMRG method is intrinsically unable to deal 
with disordered systems because randomness would in- 
validate building up a system in a block like fashion. 

Our aim in the present work is to show that a careful 
DMRG calculation can indeed be successfully applied in 
the presence of quenched disorder. In particular, we pro- 
vide a quantitative phase diagram for the Id disordered 
Bose-Hubbard model at both integer and half filling. For 



integer filling it is found that the superfluid and Mott in- 
sulating state are always separated by a Bose glass phase 
as suggested by Fisher et al.cl The superfluid density is 
nonmonotonic not only as a function of interaction but 
also of disorder. Thus for strong repulsion increasing 
disorder drives a transition from a Bose glass to a su- 
perfluid. For half filling, where no Mott insulator exists, 
the superfluid density is again a nonmonotonic function 
of the repulsive interaction, however disorder now always 
suppresses superfluidity as expected. The corresponding 
phase diagram is in agreement with that suggested by 
Giamarchi and SchulzEl, however we find no indication of 
a qualitative difference between the glass phase at small 
or large values of the repulsion (Anderson vs. Bose-glass). 

The BoseJiiibbard model in Id is defined by the 
HamiltonianErH 

H = - l - ^2(b\ +1 bi + h.c.) + ~ X ) + e *" 1 - 

i i i 

(i) 

Here b\ is the boson creation operator on site i of a Id lat- 
tice with L sites and = b\bi the corresponding local oc- 
cupation number with eigenvalues 0, 1,2,.. .. The kinetic 
energy is described by a hopping matrix element t > 0, 
leading to a standard tight binding band e(k) = —tcosk 
in the absence of interactions and randomness. The re- 
pulsive interaction is described by a local, positive Hub- 
bard U which increases the energy if more than one bo- 
son occupies a given site. Finally the site energies q are 
assumed to be independent random variables with zero 
average and a box distribution in the interval from —A 
to A. Throughout we work in the canonical ensemble 
with a given (dimensionless) density n — ^- of bosons. 
As usual we choose t — 1 as a unit of energy (note that 
some authors have t instead of % in the hopping or 2e^ in 
the site energies which leads to a trivial factor of two dif- 
ference with our results). Apart from the density n, this 
leaves the two dimensionless parameters U and A char- 
acterizing the interactions and disorder which completely 
specify the problem at zero temperature. In order to dis- 
tinguish the various possible phases, we calculate both 
the energy gap E g and the superfluid fraction p s . The 
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energy gap which is only nonzero in a Mott insulating 
phase, can either be evaluated directly from a numeri- 
cal calculation of the energy of the ground and first ex- 
cited state. Alternatively the gap can be obtained as the 
difference E g = p, p — p n between the chemical potent 
tial for particle (p p = En+i ~ En) or hole excitationstl 
(fi n = En — En-i)- We have employed both methods in 
order to check our results. For the superfluid fraction p s , 
we use the thermodwiamic definition proposed by Fisher, 
Barber and Jasnowa. It is based on defining p s via the 
sensitivity to a change in the boundary conditions be- 
tween periodic (pbc) and antipcriodic (apbc) ones. In 
one dimension, at a given density n = the superfluid 
fraction p s is thus given by (t = 1) 

or j 

Ps = ^-^[E a o PbC (L)-E p bc (L)] (2) 

where Eq(L) are the ground state energies for the spe- 
cific boundary condition. In the absence of interactions 
and disorder it is straightforward to show that p s = 1 
for arbitrary densities n as it should be. It is important 
to note that it is precisely a nonvanishing value of p s 
(in the limit L — > oo) which is the relevant criterion for 
superfluidity despite the fact that the one particle den- 
sity matrix (b\bo) decays to zero algebraically, i.e. only 
exhibits quasi-long range order. 

For the numerical calculations it is obviously necessary 
to limit the number of bosons which can occupy a given 
site. In order to be able to cover also small values of U, 
where many bosons tend to cluster at locally favorable 
sites, we have truncated our basis to m = 7 states for each 
site i which allows up to 6 bosons occupying the same 
place. We have checked carefully that our results do not 
depend on m, which was the case at least down to U ~ 
0.5. In the DMRG calculation we studied system lengths 
up to L = 50 and included up to M — 190 states. For the 
truncation error, which is one minus the density matrix 
eigenvalues X a of all M states kept in the decimation, 

M 

a=l 

we find values of the order of 10~ 10 . A very important 
point which turns out to be absolutely crucial for apply- 
ing the DMRG to disordered-systems is to apply the finite 
size ( "sweeping" ) algorithm!! After the system has been 
grown to its full length, renormalization group transfor- 
mations have not yet been able to take into account the 
full structure of disorder while working on shorter sys- 
tems. The finite size algorithm then works on the com- 
plete system, and improves results essentially in a vari- 
ational fashion. We find good convergence of both the 
gap and the superfluid fraction after several sweeps. The 
dependence on the number of states kept was compar- 
atively weak (also compared to the scattering of results 
in various realisations of disorder) such that we preferred 
to invest computational resources rather in sweeps. The 



antiperiodic boundary condition has been implemented 
by replacing the hopping energy t at one of the bonds by 
—t thus enforcing a localized twist in the phase by it. 
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FIG. 1. Phase diagram for commensurate filling n — 1. Er- 
ror bars are mainly due to the dependence of results on the 
realisations of disorder. Above a disorder strength A max — 4 
it is always energetically advantageous to destroy superfluid- 
ity in favor of a glass phase. 

For the discussion of our results we first concentrate on 
a commensurable>density n = 1, where a Mott insulating 
phase is expectedEI at sufficiently large U. In the limit of 
vanishing disorder A = the system is superfluid at small 
values of U with a superfluid fraction p s which monoton- 
ically decreases from one at U — to zero at U = U c . 
Since the transition to the Mott insulator is driven by 
phase fluctuations at apgiven density, it is a Kosterlitz- 
Thouless like transition!] very similar to the one present 
in a chain of Josephson junctions with a local charging 
energy!! Our numerical result for the critical value of U 
is U C (A — 0) = 1.92 ± 0.04 whicluis surprisingly close 
to that found in mean field theory^. It also agrees with 
a very recent DMRG calculation of the Bose-Hubbard 
model without disorder by Kiihner und Monienliil. They 
have used the condition that the exponent K character- 
izing the decay of the off-diagonal density matrix 

(bfbo) ~ V\- K/2 (4) 

in the suprj-fiuid phase takes on the value K c = 1/2 at the 
transition!^. Note that K scales like y/Uji at least in a 
Josephson junction array description which is equivalent 
to the Bose-Hubbard model at large integer densities. 
At finite disorder the Mott insulating phase is suppressed 
because the energy gap is reduced. For vanishing hop- 
ping, i.e. U —> oo effectively, the reductiona is just 2A. 
Thus in the limit of large U the Mott-insulator disappears 
of A > U/2. This is in fact the asymptotic behaviour of 
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the transition line shown in Fig. 1. For nonzero t, i.e. 
finite U the transition appears earlier, until the Mott in- 
sulator completely disappears at U < U C (A = 0) = 1.92. 
Outside the Mott-insulating phase the gap vanishes, how- 
ever at finite disorder the system need not be supcrfluid. 
Indeed calculating the superfluid fraction p s , we find that 
p s is nonvanishing only in the superfluid regime in Fig. 
1, which bends down to A = both near U = and 
U = U C (A = 0). As a consequence, at finite disorder, 
there is no direct transition from a Mott insulator to 
a superfluid|-|in agreement with the arspmcnts given by 
Fisher et al.el and Freericks and Momenta The complete 
phase diagram is shown in Fig. 1. It agrees well with 
that found by Prokof'ev and SvistunovQ using a rather 
different method apd also with the qualitative picture put 
forward by HerbutHJ. By contrast, there are strong, even 
qualitative differences with the phase diagram found by 
Pai et alB. Their failure to see the intervening Bose- 
glass between the superfluid and the Mott-insulator is 
probably related to the fact that without the sweeping 
algorithm the treatment of a disordered problem by the 
DMRG is not reliable. Regarding the general structure of 
the phase diagram shown in Fig. 1, we expect that it will 
not be qualitatively different for the two dimensional case 
(though the corresponding path integral Mpnte Carlo cal- 
culations of-i^rauth, Trivedi and Ceperleyo and also more 
recent onest3 failed to see the intermediate Bose-glass be- 
tween the superfluid and the Mott-insulator). Assuming 
that the phase diagram of Fig. 1 is indeed generic for the 
disordered Bose-Hubbard model at commensurate densi- 
ties, one can draw two general conclusions: 

(i) Since the superfluid fraction is a nonmonotonic 
function of U for a given disorder, repulsive interactions 
have a delocalizing tendency at small U but enhance lo- 
calization at large U . This is in fact a general prop- 
erty, valid also at incommensurate .densities as verified 
by Scalettar, Batrouni and Zimanyitl for n = 0.625 and 
our own results at n = 0.5. 

(ii) More surprisingly, for fixed repulsion in the range 
U > U C (A = 0) but not too large, increasing disorder 
drives a Bose-glass to superfluid transition. Thus in- 
creasing disorder may in fact favour superfluidity (see 
dash-dotted line in Fig. 1). The associated superfluid 
fraction is finite only for A > A_({7). It first increases 
with A but eventually decreases to zero again at the up- 
per boundary A + (U) This effect may be understood by 
observing that with increasing distance from the Mott 
insulator the density of mobile particle-hole excitations 
increases, thus enhancing p s . At larger values of A the 
disorder induced localization takes over and p s goes to 
zero again at the upper phase boundary A + (U). 
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FIG. 2. Phase diagram for incommensurate filling n = 0.5. 
Errors due to the dependence of results on the realisations 
of disorder are now much stronger (with the exception of the 
point at A = 0). Above a disorder strength of 0.95 we find 
that ps — within our numerical accuracy. 

For the study of a noncommensurate density, we have 
chosen n = 0..5 where the Mott-insulating phase is com- 
pletely absento. The resulting phase diagram is shown in 
Fig. 2. It exhibits a superfluid phase in a finite regime 
U-(A) < U < U+(A) of the repulsion provided the dis- 
order is below a critical value A max k, 0.95. The er- 
ror bars in the determination of the phase boundary are 
larger than in the case n = 1 because the superfluid frac- 
tion exhibits rather strong realization dependent fluctu- 
ations. This problem becomes particularly relevant in 
the limit of small U. In fact noninteracting bosons are 
a singular limit of the disordered Bose-Hubbard model 
since particles will collapse into the single level with the 
lowest 6i, which will vary between different realizations. 
For small but finite U the ground state densities are still 
rather nonuniform. Now on the basis of that, it has been 
conjectured by Scalettar et al.H that there are two quali- 
tatively different localized states, a suggestion originally 
due to Giamarchi and SchulzEI. The two phases would 
be separated by a line A C (U) above A max which meets 
the phase boundary to the superfluid at a multicritical 
point. In order to look for signatures of this boundary at 
A > A max , we have calculated the expectation value of 
the dimensionless disorder energy per particle 



(5) 



which is finite for a localized stated. Although S be- 
comes increasingly negative as U is lowered, approaching 
the limit S = — 1 at U C 1, we have found no indica- 
tion of any abrupt changes. This suggests that there is 
no quantitative distinction between an "Anderson glass" 
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for small U and a Bose-glass for larger repulsion. Verly 
likely it is only the line U — which is singular. This 
point of view is supported further by the fact that the 
phase diagram found by Prokov'ev and SvistunoMH on the 
basis of the Giamarchi and Schulz criterion^ K c = 2/3 
for the renormalized exponent in the decay of the off- 
diagonal density matrix (4) essentially coincides with our 
results. Thus for any point on the phase boundary be- 
tween the superfluid and the Bose-glass, scaling is to- 
wards A = 0, K = 2/3 even for very small U. Finally we 
have also determined the superfluid fraction as a function 
of U, which exhibits a maximum at U rs 1 — 1.5. Un- 
like the case for n — 1 this maximum does not scale to 
larger U if A is increased. For very small A the critical 
value U C (A = + ) = 3.2(2) beyond which p s vanishes in 
the presence of even a very small randomness has been 
determined by calculating the exponent K in the pure 
system and using the criterion K c = 2/3. Quite gener- 
ally, however, the numerical calculation of p s becomes 
rather difficult for small disorder. This is probably re- 
lated to the strong divergence £ ~ (-j-) 1 ^ 3 2 I K ^ Q f the 
localization length in the limit A — > near the critical 
point K c = 2/3, which follows from the integration of 
the Giamarchi and Schulz flow equations. For vanishing 
disorder p s is finite for arbitrary values of U, approach- 
ing p s = 2/tt as U — > oo where the Bose- Hubbard model 
at n = ^ is equivalent to the exactly soluble quantum 
XY-modelEj in zero magnetic field. Since K — oo in 
this limit, the localisation length in the XY- model with 

1/3 

a random local field is expected to diverge like (i) 

In conclusion we have demonstrated that the DMRG 
method can be successfully applied to systems with 
quenched disorder. The phase diagram of the Id Bose- 
Hubbard model has been determined both at integer and 
at half filling. It exhibits significant differences with ear- 
lier DMRG resultsD but essentially agrees with a very re- 
cent quantum Monte Carlo calculation!! Our conclusions 
quantitatively support the general picture for the disor- 
dered Bosc-Hubbard model developed by Giamarchi and 
SchulzEl and by Fisher et al.El. The model studied here is 
probably the simplest example for the interplay between 
interactions and disorder and as such is clearly of interest 
in itself. Experimental realisations e.g. in terms of vqx=. 
tices in a Id array of Josephson junctions with disorderLD 
or the recent suggestion that Bose-Hubbard physics may 
be relevant for cold atoms in optical latticesta, will cer- 
tainly further the interest in this model. 
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